function [ correctedRatioIm rgbMpMap data plotData scoreMatrix]...
    = runMacPigAlgorithm(imageSet, imageFeatures, correctionMethod,...
    centerChoice)
%UNTITLED4 Summary of this function goes here
%   Inputs: 
%   Raw Images after computeDirAvg and imageFeatures
%   ImageFeatures: fovRow, fovCol, macRadius, macBoundary, macRegion,
%   nonMacRegion
%   Correction Method
%   Normalize or Not
%
%   Outputs:
%   Ratio image
%   TV
%   RSCAN Data
%   RGB MP Map
%   Mag Correction
%   alpha, beta, 
%   NOTE: make the normOption feature clearer.

% Normalize if MP selected centers. Do not normalize if user selected
% centers, and instead feed a normalizedImageSet into the imageSet input.

% Correct the images based on the method chosen. Make sure to define all
% outputs possible for the excel sheet. Make output place holders if the
% method needs to. For instance, Uniform correction does not produce an
% alpha, so we need a place holder for it.
if strcmp(correctionMethod, 'Variable')
    [correctedImSet alpha beta scoreMatrix negValues magCorrection]...
        = computeRedImageBgCorrection(imageSet, imageFeatures);
    bestBlueFrac = 'N/A';
    bestOtherFrac = 'N/A';
elseif strcmp(correctionMethod, 'Uniform')
    loPercentile = 0.5;
    % Need to input un-normalized images into this
    [correctedImSet scoreMatrix magCorrection bestBlueFrac bestOtherFrac]...
        = computeUniformBgCorrection(imageSet, imageFeatures, loPercentile);
    % Need to blur and normalize the corrected set for image consistancy
    % among method choices. 
    alpha = 'N/A';
    beta = 'N/A';
    negValues = 'N/A';
elseif strcmp(correctionMethod, 'Uncorrected');
    correctedImSet = imageSet;
    magCorrection.blueCenter = 0;
    magCorrection.otherCenter = 0;
    magCorrection.blueBoundary = 0;
    magCorrection.otherBoundary = 0;
    alpha = 'N/A';
    beta = 'N/A';
    negValues = 'N/A';
    bestScore = 'N/A';
    bestBlueFrac = 'N/A';
    bestOtherFrac = 'N/A';
else 
    error('Please use a valid correction string. See help menu.')
end

% Correct for axis blips (these are when the line jumps quickly to infinite
% in the plot) by setting the region outside the macular region
% equal to one number
correctedImSet.blueIm(imageFeatures.iNonMacRegion) = ...
    mean(correctedImSet.blueIm(imageFeatures.iMacBoundary));
correctedImSet.otherIm(imageFeatures.iNonMacRegion) = ...
    mean(correctedImSet.otherIm(imageFeatures.iMacBoundary));
% Compute the ratio
correctedRatioIm = correctedImSet.blueIm./correctedImSet.otherIm;
% Correct for axis blips again by setting the region outside the macular
% region equal to one number
correctedRatioIm(imageFeatures.iNonMacRegion) = ...
    mean(correctedRatioIm(imageFeatures.iMacBoundary));

% Compute the total variation of the macular pigment map, which is the
% corrected ratio image
tvOfRatioIm = compute_total_variation(correctedRatioIm);

% Compute the radial average of the pixels in the image to plot
plotData = rscan(correctedRatioIm, 'xavg', imageFeatures.fovCol, 'yavg',...
    imageFeatures.fovRow, 'rlim', imageFeatures.macRadius, 'dispflag', 0);
initialPoint = plotData(1);
halfDegreePoint = plotData(floor(length(plotData)/8));
oneDegreePoint = plotData(floor(length(plotData)/4));
finalPoint = plotData(end);

% Indicate the red dot on the macular pigment map. Provide a check against
% rgb images
[indexMpMap, map] = gray2ind(correctedRatioIm);
rgbMpMap = ind2rgb(indexMpMap, map);

% Generate a block of pixels to highlight
FOV_BOX_SIZE = 3; %pixel box that is 6 by 6
rowRange = imageFeatures.fovRow - FOV_BOX_SIZE: imageFeatures.fovRow + FOV_BOX_SIZE; %Box width
colRange = imageFeatures.fovCol - FOV_BOX_SIZE: imageFeatures.fovCol + FOV_BOX_SIZE; %Box height

% Highlight pixels in red
rgbMpMap(rowRange, colRange, 1) = 1;
rgbMpMap(rowRange, colRange, 2) = 0;
rgbMpMap(rowRange, colRange, 3) = 0;

% Create a data array containing all valuable quantitative information. 
% This will be an output that will need to be mixed with another array of 
% qualitative information containing Patient Identifier, Method, and Red
% RGB or IR. 

% Generate a cell array of general ID tags that have their assigned value
% underneath them. This is the general quantitative information that will 
% be stored
dataTags =...
    {'TV', 'Initial Point', 'Half Degree Point', 'One Degree Point',...
    'Final Point', 'Blue Correction at Center',...
    'Other Correction at Center', 'Blue Correction at Boundary',...
    'Other Correction at Boundary', 'Alpha', 'Beta',...
    'Best Blue Frac', 'Best Other Frac'...
    'fovRow', 'fovCol'};
indexPosition = 1;
% Use a for loop to dynamically name the tags according to the method and
% center used in analysis
for string = dataTags
    charString = char(string);
    tempID = {[centerChoice ' ' correctionMethod ' ' charString]};
    dataTags(1, indexPosition) = tempID;
    indexPosition = indexPosition + 1;
end

dataQuantity = {tvOfRatioIm, initialPoint, halfDegreePoint, oneDegreePoint,...
    finalPoint, magCorrection.blueCenter,...
    magCorrection.otherCenter, magCorrection.blueBoundary,...
    magCorrection.otherBoundary, alpha, beta,...
    bestBlueFrac, bestOtherFrac,...
    imageFeatures.fovRow, imageFeatures.fovCol};

data = [dataTags;dataQuantity];

end